Heterogeneous Nucleation in the Low Barrier Regime 
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In simulations of the 2D Ising model, we examine heterogeneous nucleation induced by a small 
impurity consisting of a line of I fixed spins. As / increases, we identify a limit of stability beyond 
which the metastable phase is not defined. Approaching the limit of stability, we evaluate G(n), the 
free energy when the system contains an n-site cluster of the stable phase attached to the impurity. 
We show that the conventional definition of the nucleation barrier as the difference between the 
minimum and the maximum (at n = n*) of G(n) gives poor results when used to predict the 
nucleation time (from classical nucleation theory) and the size of the critical cluster (from the 
nucleation theorem). However, if the barrier is defined relative to a reference state that considers 
the entire configuration space of the metastable phase (i.e. all n < n*), then the predictions of theory 
are excellent. Under the latter definition we show that, contrary to expectations, the nucleation 
barrier does not vanish at the limit of stability. 



Classical nucleation theory [ll-Q (CNT) continues to 
be of fundamental importance for understanding phase 
transformations in a great variety of systems. The as- 
sumptions upon which CNT is based nominally restrict 
the theory's ability to predict nucleation rates to systems 
that are only mildly metastable and for which the nucle- 
ation barrier is large relative to kT, where T is the tem- 
perature and k is Boltzmann's constant. However, many 
interesting phase changes occur in the deeply metastable 
regime where the system is approaching a limit of stabil- 
ity and the free energy barrier to nucleation is expected 
to disappear [3, [B| . Some recent simulation studies (fHU 
find that the general features of CNT remain surpris- 
ingly robust in this deeply metastable regime, but others 
suggest that CNT breaks down [1, [t| [lj| ■ Understand- 
ing whether CNT remains applicable, or how the theory 
must be adapted, when the nucleation barrier becomes 
low remains an open question. 

The presence of a heterogeneous interface in a 
metastable system can dramatically lower the nucleation 
barrier in phase transformations such as vapor conden- 
sation and crystallization [ll], [l2j]. Consequently, het- 
erogeneous nucleation plays an important role in a vari- 
ety of phenomena including atmospheric physics Tj^, 



the use of templates to form complex structures 
171 ] and protein crystallization (18(. The basic princi- 
ples of heterogeneous nucleation involving macroscopic, 
bulk surfaces are relatively well established. However, 
in many cases the heterogeneities are microscopic in size 
and there is considerable interest in understanding how 
particle size influences the nucleation mechanism and 
Til- 22j|. especially as the barrier ap- 
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proaches kT. 

In this Letter, we study heterogeneous nucleation in 
the two-dimensional (2D) Ising model to explore the na- 
ture of the nucleation barrier on approach to the limit of 



stability of a metastable phase. We seek to clarify the def- 
inition of the barrier in this limit, and to test the degree 
to which theories (in particular CNT, and also the nucle- 
ation theorem) are able to predict the behaviour observed 
directly in this regime. The Ising system we examine was 
studied previously by Sear [2(|, who demonstrated that 
a small cluster of fixed "impurity" spins increased the 
nucleation rate significantly relative to the homogeneous 
nucleation rate. In the present study, we exploit the fact 
that by increasing the size of the impurity we can sys- 
tematically lower the nucleation barrier and also bring 
the system to a limit of stability. At the same time, this 
simple model permits a rigorous evaluation of the free 
energy barrier, critical cluster size, and nucleation rate, 
facilitating the comparison of theory with direct obser- 
vations. 

Our results are based on Monte Carlo (MC) simula- 
tions of a 2D Ising model of a ferromagnet on a L x L 
square lattice with L = 45, the same system size studied 
in Ref. [201 ] . The energy of the system in spin configura- 
tion c is given by, 
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where Sj = ±1 is the spin value of site i, J > quantifies 
the ferromagnetic exchange interaction, H is the value of 
the external magnetic field, and N — L 2 is the number 
of sites in the lattice. The sum in the first term is taken 
over all nearest-neighbor pairs of spins. We explore the 
configuration space of the system using a standard MC 
dynamics, in which one Monte Carlo step (MCS) corre- 
sponds to N spin-flip attempts, where spins are chosen 
at random. 

To induce heterogeneous nucleation, our system con- 
tains an impurity consisting of a line of I spins fixed 
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FIG. 1: Example configurations of the system in the 
metastable phase with j3J = 0.65, j3H = +0.05, and / = 7 
(a) when n — no = 21, and (b) when n — n* — 174. Impurity 
sites are shown in red. Up-spins (sj = +1) are shown in black, 
down-spins (sj = — 1) in white. 



to Si = +1. In each of our runs, we initialize all non- 
impurity spins to Si = — 1, and equilibrate the system in 
the spin-down phase at j3 J = 0.65 (i.e. 0.678 of the criti- 
cal temperature) and (3H = —0.05, where f3 — 1/kT. We 
then create a metastable state by instantaneously chang- 
ing the sign of the magnetic field, so that f3H — +0.05. 
Under these conditions the spin-down phase is metastable 
(see Fig. [1]), and the system persists in this phase until 
nucleation of the stable spin-up phase occurs. 

In the absence of an impurity, the transformation 
from the metastable to the stable phase occurs via ho- 
mogeneous nucleation, for which the free energy bar- 
rier is 27 'kT [23 j_ 24 1 and the system nucleation time is 
1.5 x 10 9 MCS[2l]7 In the presence of the impurity, the 
phase transition occurs via heterogeneous nucleation. As 
we will show below, for I > 3, the height of the heteroge- 
neous nucleation barrier is less than lAkT, and the het- 
erogeneous nucleation rate is more than 350 times faster 
than homogeneous nucleation. Homogenous nucleation 
events (i.e. events that do not involve the impurity) are 
thus rare and can be neglected in our analysis. 

To find the free energy barrier for heterogeneous nucle- 
ation, we seek to evaluate the minimum reversible work 
of formation of a critical cluster of the stable phase. 
Since homogenous nucleation can be neglected, the crit- 
ical cluster is necessarily a cluster of up-spins (i.e. sites 
with Si — +1) attached to the impurity. In the following, 
we define the "impurity cluster" as the contiguous clus- 
ter of up-spins that contains the impurity spins; thus the 
number of spins n in the impurity cluster includes the 
impurity spins themselves. Under this definition there 
can only be one impurity cluster, and so n is a system 
property (and hence an order parameter) with respect to 
which the nucleation free energy barrier may be defined. 

To define the system free energy for fixed (N, H,T,l), 
we first write the partition function of the system Z = 
^2n=i Z( n ) as a sum over the conditional partition func- 
tion Z(n) = Y] c ( n \ exp(— f3E c ). The sum in Z(n) is over 
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FIG. 2: (a) G(n) - G(l) as a function of n for I = 3 to 12. 
(b) G(n) — G m versus n for I = 3 to 11. For all curves, the 
statistical error is less than 0.02fcT. 



all system configurations c in which the impurity cluster 
consists of exactly n spins. The corresponding condi- 
tional free energy is G(n) = — kT log Z(n), which is the 
free energy of the system when it contains an impurity 
cluster of size n. 

To compute G(n), we note that the probability to ob- 
serve an impurity cluster of n spins is P{n) = Z(n)/Z. 
Consequently, the work of formation of an n-spin impu- 
rity cluster, starting from a "bare" impurity (i.e. n = I), 
is given by the free energy difference, 



G(n) - G(l) = -fcTlog^TT = -fc^log 



Z{1) 



P(l) 



(2) 



To determine G(n) in the region of the nucleation barrier, 
we impose a constraint on our MC sampling such that 
n < n max — 300. This choice of n max restricts our sim- 
ulations to the metastable phase, and to configurations 
in the vicinity of transition states to the stable phase, 
without allowing the system to explore the stable phase 
itself. 

We evaluate the equilibrium ratio P(n)/P(l) from 
our simulations, and plot the result for G(n) — G(l) in 
Fig. [D(a), for various I. For 3 < I < 11, each curve 
exhibits a maximum at n = n*, indicating the size of 
the critical cluster, and demarcating the boundary be- 
tween the configuration space of the metastable phase 
(defined by I < n < n*) and the stable phase (n > n*). 
Also, as I increases, a minimum at n = no emerges and 
grows; this feature corresponds to wetting of the impu- 
rity by a finite cluster of the stable phase. The variation 
of n* and no with I is shown in Fig. [3] Note that in 
Fig. EKa), for I = 12, the variation of G(n) - G(l) is 
monotonic in n, demonstrating that the limit of stability 
of the metastable phase has been exceeded. This limit of 
stability is also seen in Fig. [3] as n* and no approach one 
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FIG. 3: n* and no as a function of I. Statistical errors are 
smaller than the symbol size. 



another, and then cease to be defined for I > 12. 

To obtain the free energy barrier for nucleation from 
G(n), we must identify an appropriate reference state 
with respect to which the barrier height will be measured. 
When both a minimum and a maximum appear in a free 
energy profile, previous works (see e.g. Refs. 0, [l2, 17|) 
have defined the free energy barrier for nucleation as 
AGo = G(n*) — G(n ), i.e. the difference in G(n) from 
the minimum at n = no to the maximum at n — n* . Our 
results for AGo as a function of / are shown in Fig. @|a). 
Under this definition, the nucleation barrier goes to zero 
at the same I that corresponds to the loss of stability of 
the metastable phase. 

However, the free energy barrier for nucleation is more 
accurately defined as the minimum reversible work re- 
quired to apply a constraint that confines the system to 
the transition state at n = n* , starting from a refer- 
ence state that considers the entire configuration space of 
the metastable phase, i.e. all configurations in the range 
I < n < n* . To implement this definition in the present 
system, we define the partition function of the metastable 

phase Z m = J2n=i ^( n ) as a restricted sum over all states 
such that I < n < n*. The corresponding free energy of 
the metastable phase is G m = — kT log Z m . The work 
of formation of an n-spin impurity cluster, starting from 
the equilibrium metastable phase, is then given by the 
free energy difference, 



G{n)-G m = -kT log 
= -fcTlog 



Z(n) 



P(n) 



(3) 



The second equality above emphasizes that G(n) — G m 
can also be evaluated in our simulations from the relative 
probabilities for observing the impurity cluster to have 
various n. 

Our results for G(n) — G m are plotted in Fig. HJb). 
The difference between the free energy curves in Fig.^a) 
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FIG. 4: (a) Comparison of AG and AGo as a function of 
I. (b) tcntj tcnt, tmfpt, and rscar as a function of I. For 
comparison, note that the homogeneous nucleation time when 
no impurity is present is 1.5 x 10 9 MCS [2(J. For all quantities, 
statistical errors are smaller than the symbol size. 



and (b) is just a change in the reference state, and so 
the shapes of the curves do not change. However, the 
height of the nucleation barrier does change. Using the 
more accurate definition, the work of formation of the 
transition state from the metastable phase is given by 
AG = G(n*) - G m . As shown in Fig.H^a), AG does not 
go to zero at the limit of stability. Although paradoxical 
at first glance, this result is physically reasonable for our 
system. Since n* remains non-zero even at the limit of 
stability, the metastable phase encompasses a consider- 
able region of configuration space (I < n < n*) up to the 
point where stability is lost. Hence the work required 
to create the transition state remains finite, even as the 
metastable state ceases to exist as a distinct phase. In 
previous work, it has been assumed that the nucleation 
barrier should go to zero as the thermodynamic stability 
of a metastable phase is lost 0, HI • Our system provides 
a counter-example. 

We next assess the implications of our results when 
estimating the nucleation time using CNT 24| . For our 
system, the CNT prediction for the nucleation time is, 



TCNT 



(/+z)- 1 exp(/?AG), 



(4) 



where tcnt is the average time (in MCS) per impurity 
for a critical cluster to appear in the system that sub- 
sequently evolves into the stable phase, z = y / f3rj/2iT 
is the Zeldovich factor, where ij = —(d 2 G/dn 2 ) n=n * is 
the curvature of G(n) at the top of the barrier. We esti- 
mate 77 from a quadratic fit to data that lies within 0.2kT 
of the maximum of G{n). fj is the attachment rate of 
monomers to the critical cluster. We determine from 
the time dependence of fluctuations of the size of critical 
clusters, following the same procedure used in Ref. [2~3 |. 
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FIG. 5: Comparison of AnNT, An NT , and An s i m 
as a function of /3H, for the I = 7 system. 

The result for tcnt obtained from our data is shown in 
Fig.BJb); we also show Tq NT , the CNT prediction for the 
nucleation time obtained if we use AGo instead of AG in 
Eq. |U We find that for the lowest barriers (at large I), 
i~ cnt and Tq NT differ by more that two orders of magni- 
tude, emphasizing the importance of selecting the correct 
definition of the free energy barrier in this regime. 

To test the accuracy of tcnt, we directly evaluate the 
nucleation time in terms of the mean first passage time 
(MFPT) for the impurity cluster to grow to the critical 
size. For a given I, we set n max = n* so that the sys- 
tem is confined to explore only the configuration space of 
the metastable phase, and bring this constrained system 
into equilibrium. Then, at a randomly selected time, we 
set t = and measure the time it takes for the system 
to first reach n = n max . The MFPT is the average of 
many such measurements. We define the nucleation time 
tmfpt as twice the MFPT, because only half of the runs 
that reach the transition state would ultimately evolve 
into the stable phase. As shown in Fig. 0Jb), tmfpt is 
in excellent agreement with tcnt, but is clearly distinct 
from Tq NT . Fig. HJb) also shows that tmfpt is consistent 
with the nucleation times (rs ea r) reported in Ref. [20] for 
the same system, as found using a forward-flux sampling 
method. Our observations not only validate the choice of 
AG over AGo for the nucleation barrier, but also demon- 
strate that in our case CNT is capable of predicting the 
nucleation time with remarkable accuracy even at the 
very limit of stability of the metastable phase. 

Finally, it is also possible to test the validity of AG 
versus AGo as a measure of the nucleation barrier by 
exploiting the nucleation theorem. The nucleation theo- 
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= -An, 



(5) 



number of up-spins in the critical cluster. In the first 
equality, we have used Afj, « 2H, which for the Ising 
model is a good approximation for T below the Curie 
temperature 13, HI- To conduct this test, we carry out 
new runs for the case of I — 7 over a range of j3H from 
0.048 to 0.056. We directly evaluate An as the difference 
in the average number of up-spins in the entire system 
(including those not in the impurity cluster) when the 
system is at n = n* , and the average number of up-spins 
in the metastable phase averaged over all I < n < ri*; 
these results are shown in Fig. [5] and denoted as An S j m . 
We also evaluate AG and AGo as a function of H, and 
estimate the derivative in Eq. [5] using a five-point central 
difference numerical method. We thus obtain two esti- 
mates of An from the left-hand side of Eq. [5) AnNT is 
found using AG, while An^ T is found using AGo- As 
shown in Fig. [5j AnNT is in distinctly better agreement 
with An s i m than An^j T , providing additional confirma- 
tion that AG is the correct measure of the nucleation 
barrier. 

We note that although the definitions of G(n) and the 
nucleation barrier are straightforward in the present sys- 
tem, it remains a challenge to define the free energy of 
cluster formation in simulations of homogeneous nucle- 
ation when the barrier is low. In particular, in the ho- 
mogeneous case, the cluster size n can only be used as 
an order parameter in the limit that the spontaneous for- 
mation of clusters is rare 



32. 33 , i.e. 



where A/i is the difference in chemical potential between 
the stable and metastable phases, and An is the excess 



in the limit that 
there is at most one cluster in the system. However, 
Z(n) [and hence G(n)] is not defined if n does not have 
a unique value for each system configuration, a compli- 
cation that arises when the barrier is low because mul- 
tiple clusters may form simultaneously. This difficulty is 
avoided here because homogeneous nucleation is negligi- 
ble, and because a single impurity cluster is always well 
defined, regardless of the barrier height. 

We also emphasize that the difference between AG 
and AGo is only apparent in the low barrier regime. 
When the barrier is high, even small clusters are rare, 
and the properties of the metastable phase are domi- 
nated by system configurations found near n = I. In 
this limit G(l) « G m , and AG and AGo become equiva- 
lent. However, as we have seen, in the low barrier regime 
the distinction between AGo and AG becomes impor- 
tant. Our analysis shows that the selection of the correct 
thermodynamic reference state is crucial for identifying 
the nucleation barrier (AG) that is appropriate for study- 
ing low barrier systems. Using AG, we demonstrate that 
both CNT and the nucleation theorem are impressively 
accurate, even at the limit of stability. At the same time, 
we find that AG does not vanish at the limit of stabil- 
ity of the metastable phase. We expect that the pattern 
of behavior found here will be common to all low-barrier 
systems where a free energy minimum and maximum con- 
verge at a finite value of the order parameter, and thus 
may be generic for heterogeneous nucleation on small im- 
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purities. Only in the case that the size of the critical 
cluster goes to zero (n* — > 0) should we expect that the 
nuclcation barrier will vanish at the limit of stability. 
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